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Condensation of a liquid droplet from a supersaturated vapour phase is initiated by a prototypical nucleation 
event. As such it is challenging to compute its rate from atomistic molecular dynamics simulations. In fact at 
realistic supersaturation conditions condensation occurs on time scales that far exceed what can be reached 
with conventional molecular dynamics methods. Another known problem in this context is the distortion of 
the free energy profile associated to nucleation due to the small, finite size of typical simulation boxes. In 
this work the problem of time scale is addressed with a recently developed enhanced sampling method while 
contextually correcting for finite size effects. We demonstrate our approach by studying the condensation of 
argon, and showing that characteristic nucleation times of the order of magnitude of hours can be reliably 
calculated, approaching realistic supersaturation conditions, thus bridging the gap between what standard 
molecular dynamics simulations can do and real physical systems. 


I. INTRODUCTION 

Nucleation is the event initiating first order phase tran¬ 
sitions in which a small embryo of a thermodynamically 
stable phase appears within a parent metastable phase. 
The formation of gas bubbles in a liquid, liquid droplets 
in a vapour or crystal particles in solution are all exam¬ 
ples of nucleation events playing a key roles in a variety 
of fields ranging from atmospheric physics to pharmaceu¬ 
tical manufacturing. The small length scale characteris¬ 
ing nucleation events renders their direct experimental 
observation inherently challenging, while providing the 
ideal playground to apply and develop molecular mod¬ 
elling techniques. Despite extensive efforts in the inves¬ 
tigation of nucleation phenomena with molecular simu¬ 
lations, the development of a systematic approach to the 
calculation of nucleation rates from first principles still 
remains a challenge, due to the very nature of the nucle¬ 
ation phenomena. 

In the context of the condensation of a liquid phase 
from a supersaturated vapour, the nucleation rate J is 
defined to be the number of liquid droplets formed per 
unit time and volume. 

The formation of a liquid droplet containing n 
molecules, in a system at constant volume (V) and tem¬ 
perature (T), is associated to a Helmholtz Free Energy 
change A F(n). The maximum of this quantity, AF*, 
corresponding to a critical number of molecules nj, con¬ 
stitutes the energy barrier that the system has to over¬ 
come in order to undergo the nucleation process. The 
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very existence of this barrier determines the key features 
of nucleation, namely that it is an activated process and 
a paradigmatic example of a rare event. 

The free energy barrier AT* and the nucleation rate 
J depend on the thermodynamic driving force, that is 
typically expressed in terms of the supersaturation S, i.e. 
the ratio of the actual vapour pressure and the equilib¬ 
rium vapour pressure. They depend also on temperature, 
T and on system-specific properties, namely the surface 
tension between the liquid and vapour phases, the molec¬ 
ular volumes of the two phases, V£ and v g , and the spe¬ 
cific surface of the newly formed droplets, i.e. 6/d for a 
spherical droplet of diameter d. 

In a system of volume V, at supersaturation S, and 
temperature T the characteristic time of nucleation r can 
be expressed as: 


J(5, T) V v ' 

In Eq. 1 it can be readily seen that, at any given con¬ 
dition of temperature and supersaturation, r, i.e. the 
average time necessary to observe a nucleation event, 
is inversely proportional to V. This relationship repre¬ 
sents a constraint between two key factors determining 
the computational cost of MD simulations of nucleation 
events: the number of steps required to observe a nucle¬ 
ation event that needs to be of order of r/St where St 
is the integration timestep used in MD, and the number 
of degrees of freedom which is instead proportional to 
the number of atoms, and thus to the volume V when 
comparing systems at the same density p = N/V. Such 
a constraint limits the range of conditions that could be 
directly investigated by small-scale MD simulations to 
regimes where J ^ 10 22 — 10 25 cm -3 s _1 . Up to now, 
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overcoming such limitations has only been possible by 
using large scale simulations involving millions of atoms 
and requiring massive computational resources 1 . 

Nucleation is the prototypical example of a rare event, 
inherently stochastic in nature. When explicitly consid¬ 
ering the stochastic character of nucleation r represents 
the expected time for the nucleation of the first liquid 
droplet from a supersaturated vapour. The law of rare 
events suggests that the first nucleation event can be in¬ 
terpreted as a Poisson process 2 ,where the survival proba¬ 
bility Po, i.e. the probability that at a given time t there 
are no droplets in volume V , is: 

P 0 (t) = exp ( 2 ) 

Eq. 2 highlights how in order to reliably estimate r, 
the stochastic nature of the nucleation process needs to 
be explicitly considered, while the reconstruction of the 
distribution of transition times is pivotal to reliably esti¬ 
mating J. In addition to the timescale issues, molecular 
simulations of nucleation processes suffer from intrinsic 
finite-size effects that cause a systematic distortion of 
A F(n) and impact rate calculations 3-8 . 

In this work we propose a systematic approach for the 
calculation of nucleation rates from small-scale unseeded 
molecular simulations. The proposed approach is based 
on recent developments of Well Tempered Metadynamics 
(WTmetaD) 9-11 that enable for the calculation of tran¬ 
sition times distributions from biased simulations. More¬ 
over we develop a systematic correction for the effect of 
finite-size on rate calculations. 

Our method is tested on the paradigmatic case of the 
nucleation of a liquid argon droplet from a supersatu¬ 
rated argon vapour. The choice of this system has a 
two-fold aim: the first is to analyse a simple and yet very 
significant system that allows to draw conclusions gen¬ 
eral enough to be extended to a wider class of nucleation 
phenomena. The second is to point out that even in such 
a simple system, at realistic values of supersaturation the 
nucleation timescale rapidly grows out of reach of stan¬ 
dard MD. Moreover choosing argon as a test case allows 
to benchmark our results against the existing literature 
based on small-scale unbiased MD simulations 12 . 

The paper is structured as follows; at first the details 
of the application of WTmetaD to the calculation of nu¬ 
cleation rates are reported, then an analysis of finite-size 
effects is carried out. WTmetaD and the finite size cor¬ 
rection are then combined to outline a systematic strat¬ 
egy for the calculation of nucleation rates. Finally results 
are reported and commented. Unless otherwise noted, 
the subscript N will be used to refer to relevant quanti¬ 
ties in finite-size systems. Such subscript will be dropped 
whenever referring to their counterpart in macroscopic 
systems. 


II. NUCLEATION RATES AND LONG TIMESCALES 
A. From metadynamics to dynamics 

In this work the acceleration effect associated to WT¬ 
metaD has been exploited in order to substantially re¬ 
duce the simulation time required to observe a nucleation 
event while simultaneously maintaining the system size 
small, hence significantly diminishing the overall compu¬ 
tational cost. WTmetaD is conventionally used to com¬ 
pute free energy surfaces in a variety of contexts 10,11 . 
Recently it has been shown that, taking inspiration from 
conformational flooding 13 and hyperdynamics 14 , transi¬ 
tion times associated to activated events can be efficiently 
computed from WTmetaD simulations 15 . In WTmetaD 
the simulated system evolves in a transformed time co¬ 
ordinate, twT, due to the application of the history- 
dependent bias potential Vb(£,£) constructed as a func¬ 
tion for the collective variable £ 10,15 . As discussed in 
detail in Ref. 15, rate calculations via WTmetaD do not 
require a converged estimate of free energy profiles, be¬ 
ing instead based on the systematic evaluation of the 
so-called acceleration factor , which represents the ratio 
between the physical time and the meta dynamics time. 

In the context of a nucleation problem, the specific 
transition time associated to a nucleation event t n can be 
computed from the corresponding WTmetaD simulation 
time t n wT as: 

tn = U,lTT(eX p (/3 Vb(£, t)))wT (3) 

where the term (exp (/3 Vb(£, t)))wr is the acceleration 
factor called a in the following. Note that (3 = 1/ksT. 
Crucial to this procedure is the hypothesis of negligible 
bias deposition at the transition state. To comply with 
such hypothesis, the bias potential is constructed through 
the infrequent deposition of potential Gaussians, in a 
properly defined space of collective variables. The fulfill¬ 
ment of such condition can be checked a posteriori using 
the approach detailed in Ref. Salvalaglio2014assessing. 
When this is the case, the whole transition time distri¬ 
bution can be recovered from a set of WTmetaD simu¬ 
lations. This approach has been applied to several prob¬ 
lems, thus allowing the computation of rates of activated 
processes such as DNA unfolding 16 , and protein-ligand 
unbinding 17 19 . 


B. A collective variable to describe liquid argon nucleation 

In WTmetaD the bias potential Ve(£, t ) is constructed 
as a function of a collective variable £ 10 db 20 j n work 
we choose as collective variable n, the total number of 
liquid argon atoms in the system, corresponding to the 
typical reaction coordinate used in classical nucleation 
theory (CNT). In order to be used in the framework of 
WTmetaD, n has been expressed as a continuous and dif¬ 
ferentiable function of the atomic coordinates 9,10,21 . To 



3 


this aim the ten Wolde-Frenkel definition 12,22 has been 
applied, in which atoms are considered liquid when they 
possess a coordination number larger than a threshold 
value q, that is chosen to be 5. The coordination num¬ 
ber of each molecule in the system is defined in a con¬ 
tinuous and differentiable form through the expression: 
c i = where is the Cartesian distance be¬ 

tween atoms i and j and /(r^) is the switching function: 


f( r ij ) = 


1 - ( nj/r c f 

1 - (rij/r c ) 12 


( 4 ) 


The number of molecules possessing <7 > q is thus cal¬ 
culated using the same functional form as in Eq. 4: 


N 1 

= £ — 

4-r l- 


(Cf/Cj ) 6 


(Cl/Ci) 


12 


( 5 ) 


C. Detecting nucleation events in WTmetaD simulations 

The biased transition time t n yvT is the simulation time 
associated with the occurrence of a nucleation event in 
a WTmetaD simulation. Hence to compute t n yvr , it 
is necessary to reliably detect nucleation events. In our 
case we apply the approach described in Ref. 12 which is 
based on the fact that a clear timescale separation exists 
between the residence time in the supersaturated vapour 
state and the time necessary for a supercritical nucleus to 
grow in size. The latter phenomenon is orders of magni¬ 
tude faster than the former, rendering nucleation a rare 
but fast event. As done in Ref. 12 , the nucleation time 
tn,wr can be directly calculated from the time evolu¬ 
tion of n(£), as the simulation time needed to overcome 
a threshold size of the emerging liquid droplet n. This 
approach remains valid as long as the threshold size cho¬ 
sen to define the transition criterion is larger than the 
critical nucleus. In such case the transition time can be 
safely considered independent of the specific value of n 12 . 


D. Expected nucleation time and nucleation rates 

As briefly mentioned in the introduction, due to the 
activated nature of nucleation, its transition time prob¬ 
ability distribution is described by the so-called law of 
rare events , and is thus expected to be exponential. The 
nucleation process, particularly the formation of the first 
nucleus, i.e. the event that matters at the scale of the 
MD simulation box, can in fact be modelled by a time- 
homogeneous Poisson process characterised by a survival 
probability Po(t) = exp(— t/rjsr). Due to its inherent 
stochasticity, an appropriate sampling of the nucleation 
times distribution is required to evaluate the expected 
characteristic nucleation time r^v 12,23,24 . In order to com¬ 
pute the expected nucleation time tn in a finite-size sys¬ 
tem, we perform a large number of independent WT¬ 
metaD NVT simulations and extract from each of them 


a nucleation time t^- After this, the survival probability 
distribution constructed from the tjy values is analysed 
and its statistical compatibility with a Poisson process 
quantified. This allows to check whether the conditions 
under which Eq. 3 is valid are satisfied 25 . The character¬ 
istic nucleation time fitted from the survival probability 
distribution tn is thus used to compute the nucleation 
rate in the finite size system as Jn = ('LvY) -1 , where V 
is the system’s volume 12,23,24 . 

E. WTmetaD simulation details 

Transition times were computed from NVT simulations 
of systems consisting of 512 argon atoms at a constant 
temperature of 80.7 K. A Lennard-Jones potential with 
e = 0.99797 kJ/mol and a = 0.3405 nm has been adopted 
to describe the interactions between argon atoms 12 . The 
potential has been truncated, but not shifted, with a cut¬ 
off length of 6.75cr. The time step for the integration of 
the equations of motion has been set to 5 fs 12 . These 
specific conditions were chosen in order to compare our 
results with those reported in Ref. 12. The equilibrium 
vapour pressure (p e ) of argon under these conditions is 
equal to 0.43 bar 12,26 . To carry out rate calculations in 
a wide supersaturation range we have followed the ap¬ 
proach proposed in Ref. 12, carrying out a series of NVT 
simulations in cubic boxes of increasing volume, corre¬ 
sponding to supersaturation values ranging from 11.4 to 
5.4. 


TABLE I. WTmetaD simulation setup summary. 


Label 

s 

1 

[nm] 

P 

[bar] 

n s im 

a 

wo At 

[kJ/mol] [ps] 

7 

Vg(nl) 

[Y/N] 

Si 

11.4 

10.5 

4.86 

100 

0.5/1.0 

0.01 

25 

5 

N 

s 2 

8.68 

11.5 

3.70 

100 

0.5 

0.01 

25 

5 

N 

S 3 

6.76 

12.5 

2.88 

50 

1.0 

0.02 

25 

5 

N 

S 4 

6.01 

13.0 

2.56 

50 

0.5 

0.01 

25 

5 

Y 

s 5 

5.36 

13.5 

2.28 

50 

0.5 

0.01 

25 

5 

Y 


In Tab. I setup parameters of the five sets of simula¬ 
tions, namely the supersaturation level 5, dimension of 
the simulation box edge /, initial pressure p, and num¬ 
ber of independent simulations per supersaturation level 
Usim have been reported. In Tab. I we also report the 
WTmetaD setup parameters, namely width of the de¬ 
posited Gaussians S, their initial height cjo , the deposi¬ 
tion stride At and the 7 factor. For a detailed description 
of the metadynamics algorithm and parameters the in¬ 
terested reader is invited to check the references 9-11 . In 
the last column we indicate whether an initial bias po¬ 
tential Vg(n) was applied (Y) or not (N) (see the results 
section for a description of simulations with and with¬ 
out Vg(n)). The highest supersaturation at which we 
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have performed simulations is S = 11.4, corresponding 
to the lowest supersaturation at which the standard sim¬ 
ulations of Ref . 12 were performed. This allowed us to 
check that our simulation setup was correctly reproduc¬ 
ing nucleation rates both in biased and unbiased simula¬ 
tions. For each supersaturation the survival probability 
distribution has been constructed by performing 50 to 
100 independent nucleation simulations. The bias was 
updated every 5000 integration steps, with a bias-factor 
of 5, and a Gaussian height of 1 x 10 -2 or 2 x 10 -2 
kJ/mol. The values of the collective variable n and of 
the total bias V B (n,t) have been collected every 100 
steps. Temperature has been controlled using the Bussi- 
Donadio-Parrinello thermostat 27 , with a time constant 
of 0.1 ps. WTmetaD simulations were performed with 
Gromacs 4.6.3 28 equipped with PLUMED 2.0 29 . 


III. NUCLEATION RATES AND FINITE SIZE EFFECTS 


The nucleation rate of droplets from a vapour can 
be explicitly derived within classical nucleation theory 
(CNT) as: 


J = A exp (-0A F*) ( 6 ) 

where A is a pre-exponential factor, p is the pressure in 
the vapour phase, p e the equilibrium vapour pressure, 
and A F* the free energy barrier to nucleation. 



FIG. 1. Nucleation of a liquid argon droplet from supersatu¬ 
rated vapour, top Free energy profiles predicted by CNT in an 
infinitely large system at constant supersaturation (A F(n)) 
and for a finite-size, confined system (AFjv(n)). Both free 
energy profiles refer to a system of 512 argon atoms, in a 
volume of 2197 nm 3 , at supersaturation £=6, with a surface 
energy cra=9 k B T. bottom Representation of the system in 
its vapour V and droplet D configurations. 


Equation 6 can be viewed as the product of two distinct 
contributions: an energetic part, corresponding to the ex¬ 
ponential term, and a kinetic one, related to the molec¬ 
ular collisions and given by the pre-exponential term. In 
order to derive a systematic correction to nucleation rates 
computed from small-scale molecular simulations, in the 
following we shall assess the impact of the finite size of 
the simulation box on both terms. 


A. Free energy of nucleation in a confined system 


We define confinement as the impossibility of exchang¬ 
ing atoms between the system and the surrounding en¬ 
vironment. Under this definition a NVT simulation 
box, represents a prototypical confined system, as its to¬ 
tal number of atoms N is by definition constant. As 
highlighted in several reference works 3 8 , the free energy 
change associated with a nucleation process is affected by 
confinement. In Fig. 1, the comparison between nucle¬ 
ation free energies in finite-size (blue) and macroscopic 
(red) conditions is illustrated together with the represen¬ 
tation of typical configurations of Argon atoms in the 
vapour (V) and liquid droplet (D) states. Hereafter we 
shall summarise how such a distortion in the free en¬ 
ergy profile affects nucleation rates 3 . We treat the Argon 
vapour as an ideal gas hence the pressure of the vapour 
before droplet formation (initial state of the system) is: 


Po(N,V,T) 


Nk B T 

V 


( 7 ) 


A liquid spherical embryo of n molecules has a volume 
Vd, and a surface Ad, that can be expressed as a function 
of n as: 


V d = nvi 


(8) 


A d = (an ) 2/3 


( 9 ) 


where V£ is the molecular volume in the liquid phase, and 
a = 67 r 1 / 2 ^. After formation of such an embryo (the 
final state of the system) the vapour pressure attains the 
following value: 


p(N,n,V,T) 


(N - n) k B T 
V — nvi 


( 10 ) 


The Helmholtz free energy change, A Fn, for the transi¬ 
tion from the initial to the final state, i.e. for the for¬ 
mation at constant temperature T and volume V of a 
n-molecule droplet from a vapour consisting initially of 
N molecules, is given by 4 : 


A F N (N,n,V,T) 


— n/3 1 In ^ + 7 (an ) 2 / 3 + 

+ IV /? -1 In (A^j + n (/? -1 - v t p e ) 

( 11 ) 
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FIG. 2. Free energy of nucleation in a finite-size system at constant N=512, T=80.7 K, as a function of the system volume 
V. The solid blue line represents the locus of the maxima of A Fn, corresponding to the critial nuclei. The dashed blue line 
represents the locus of the local minima in of A Fn, representing a stable argon droplet in a finite-size induced equilibrium 
with the argon vapour. In solid red the A F(n) is highlighted, corresponding to the threshold value of V above which A F(n) 
becomes a monotonically increasing function even if aSq > 1- 


where volume, surface, and pressure effects are accounted 
for. A typical AFy(n) profile is reported in Fig. 1 (a) 
(blue). Both a local maximum and a local minimum can 
be identified along AFy(n). The local maximum A F^- 
represents the free energy barrier to nucleation in a finite 
size system. Ais associated with a critical nucleus size 
n] y, that represents the liquid embryo in unstable equi¬ 
librium with the surrounding vapour. Such Avalue 
can be computed numerically. As extensively discussed in 
Ref. 3 for an argon vapour and in Refs. 7, 8 for the case of 
crystal nucleation from solution, there exists a minimum 
value of the initial supersaturation, So = po/p e > 1? be¬ 
low which the function AFy(n) is monotonically increas¬ 
ing, hence no maximum is present. At N and T fixed 
such a condition defines an upper bound for the volume 
for which nucleation rates can be computed. In Fig. 2 
AFy(n) is plotted as a function of the system volume, 
highlighting the critical conditions associated with the 
transition to a monotonically increasing function. Con¬ 
trary to the analysis so far, Classical Nucleation Theory 
(CNT) deals with infinitely large systems, where the for¬ 
mation of the liquid droplet has a negligible effect on 
the state of the surrounding phase. The corresponding 
Helmholtz free energy change, AF(n), for the formation 
of a n-molecule embryo at temperature T and supersat¬ 
uration So can be obtained by taking the limit of Eq. 
11 with V and N approaching infinity, and their ratio 
remaining constant. Under these conditions p = po and: 


AF(n) = — nkBTlnSo + 7 (an) 2 ^ 3 (12) 

For an infinitely large system, in which Eq. 12 holds 
the free energy barrier AF* can be computed analytically 


4/3 3 7 3 a 2 
27 (In Sq) 2 


(13) 


It is worth noticing that Awhich is the nucleation 
free energy barrier in a confined system at the same con¬ 
ditions of T and So, is strictly larger than A F*. 


B. Macroscopic nucleation rates from finite size 
calculations 

In order to compute a correction term associated to the 
confinement effect we follow the approach of Ref. 3, and 
define a factor as the ratio between the nucleation rate 
in macroscopic conditions J and in a finite sized confined 
system J^. In analogy with Eq. 6, the nucleation rate 
in a confined system Jn can be written as 4 

Jn = AnS exp (— /3AFtf) (14) 

Since the system supersaturation in a finite-size simula¬ 
tion before a supercritical nucleus forms is essentially the 
same as in a macroscopic system, S = So, <t> reduces to 3 : 

4> = ~r~ = ~j~ ex P (/5(Ai^ — AF*)) (15) 

J N A/v 

Eq. 15 provides a working principle to obtain nucle¬ 
ation rates in macroscopic systems from finite size NVT 
simulations as: 


J = 4>Jn (lb) 

As reported in Ref. 3, due to the exponential depen¬ 
dence on the strictly positive quantity A F^ — AF *, the 
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FIG. 3. Workflow summary for the calcuations of macroscopic nucleation rates from finite size NVT WTmetaD nucleation 
simulations. 


dominating term in Eq. 15 is exp (/3(AF^ — AF*)), 
whereas the pre-exponential term A/An in Eq. 15 is 
expected to play a secondary role. In order to identify 
key contributions to A/An, its dependence on finite size 
is discussed in the following. 

Within the framework of CNT, A is typically expressed 
as 12 ’ 30 

A = zf *^f ^ 

where Z is the Zeldovich factor given by 31 : 


as: A Fjy, A F*, Z/v, Z , n^y, and n*. Both A F(ri) and 
A F(n) N can be respectively computed from Eq. 12 and 
11 , once the surface tension 7 is known. 

The surface tension 7 is obtained by fitting Eq. 14, 
in which the pre-exponential term An in considered 
supersaturation-independent, on the Jn values obtained 
as a function of supersaturation. 

IV. WORKFLOW SUMMARY 


Z = 


I d 2 AF(n ) I 
| dn 2 \n=n* 

2irkBT 


(18) 


and /* the rate of attachment of molecules to the critical 
cluster. Since nucleation of a droplet from its vapour is a 
process controlled by direct impingement 30 , the attach¬ 
ment rate /* is derived from the kinetic theory of gases 
as 30 ’ 32 : 


/* = c(n*) 


P 

\f 7 hvmk^T 


(19) 


where c(n*) = \J(36ttV£ 2 )(n*) 2 / 3 is the surface area of 
the critical cluster, V£ is the volume per molecule in the 
liquid phase and p the pressure. 

The attachment frequencies /* and differ due to 
two reasons. The first is that the critical nucleus size 
in finite size simulations is strictly larger than the 
critical nucleus in the corresponding infinite case n* 4 . 
The second reason is that the vapour pressure acting on 
the critical nucleus in finite size systems p ^ = (N — 
n* N )/((y — n* N vi)kBT ) is always smaller than its corre¬ 
sponding value for a system at macroscopic conditions 
p = N/{Vk B T). 

The Zeldovich factors Z and Zn are instead expected 
to differ due to the fact that the curvature of the free 
energy profile in the region around its maximum is af¬ 
fected by finite size, see Fig. 1 for example. The extent 
of the contribution of the exp (/3(AF^ — AF*)), /*//am 
and Z)Zn terms to (j) is discussed in the Results section. 

The correction factor </> depends on quantities that can 
be directly calculated from A F(n) N and A F(n) such 


In Fig. 3 the workflow for the calculation of nucle¬ 
ation rates from small-scale finite size NVT simulations 
has been summarised. The calculation procedure can be 
outlined as follows: 

1. A set of WTmetaD simulations is carried out for 
multiple supersaturation levels. Supersaturation is 
imposed by defining the system volume while keep¬ 
ing constant the number of molecules N and the 
temperature T. 

2 . Applying the criterion for the identification of nu¬ 
cleation events proposed in Ref. 12, the WTmetaD 
transition time twr and the corresponding acceler¬ 
ation factor a are calculated from each WTmetaD 
simulation. 

3. The physical transition time associated to each nu¬ 
cleation event is computed using Eq. 3. 

4. The transition times obtained for each supersatu¬ 
ration value are used to fit the survival probability 
distribution, and compute the average nucleation 
time tn for each finite size system at volume V 
and supersaturation S. 

5. Average nucleation times tn are converted to finite- 
size nucleation rates using J N = 1/ (r^V) 12 . 

6 . The finite-size nucleation rates are used to fit Eq. 
14. The fitting parameter is the surface tension 7 , 
which is used to compute AF*, AF^-, n*, Z, 
Zn, and thus the correction factor <fi. 
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7. Eq. 16 is used to compute the nucleation rate in 
macroscopic conditions J . 


V. RESULTS 
A. WTmetaD simulations 

The time evolution of the number of liquid-like argon 
atoms n in a typical WTmetaD simulation is reported in 
Fig. 4, where the nucleation event can be clearly iden¬ 
tified as the rapid transition from n values fluctuating 
close to zero, to n values fluctuating around a positive 
value rip • The final state corresponds to a finite-sized 
droplet stabilised by finite size effects corresponding to 
the local minimum in free energy shown in Fig. 1 (blue 
curve). 



FIG. 4. Time series of the collective variable n{t) obtained 
from a typical simulation. The value of typical n, n* N , and 
no have been highlighted on the plot. 

It can be seen that the lifetime of the supersaturated 
vapour state in the WTmetaD simulation is much larger 
than the transition time associated with the nucleation 
event driving the system into the stable state charac¬ 
terised by n = tid • Such a difference becomes ex¬ 
ponentially large when the WTmetaD time is rescaled 
to real time according to Eq. 3 15 , neatly highlighting 
the timescale separation characteristic of the nucleation 
problem. During WTmetaD simulations a repulsive bias 
potential is adaptively constructed with an infrequent de¬ 
position of Gaussians 9,15 . In order to speed up the adap¬ 
tive construction of the bias for the two slowest cases (S 4 
and S§ in Tab . I), in addition to the WTmetaD bias 
Ve(n,t), we apply a static bias V^(n) constructed from 
a preliminary WTmetaD simulation. In Fig. 5 the total 
bias potential is reported for supersaturation levels S 2 
and $ 5 , which are characterised by the absence and pres¬ 
ence of an initial bias V#(n), respectively. In all cases 
in the region of the maximum of AjyF(n), the total bias 
applied V^ ot (n) decays to values smaller than /c^T, in 
agreement with the hypotheses of negligible bias depo¬ 
sition at the transition state invoked to carry out rate 


calculations from WTmetaD 15,25 . 



n 



FIG. 5. a) Simulation set S 2 : WTmetaD bias potential 
at transition time, Ve(n, £vct). b) Simulation set S4 ' ini¬ 
tial static bias potential Vb(ra), WTmetaD bias potential at 
transition time VB(n,twT ), and total bias at transition time 
Vi 3 0t ( n ,twT) — Vb(ti) + Vb^ti^wt)- For comparison nucle¬ 
ation free energy profiles A Fn(ji) have been reported and the 
ksT level has been highlighted in red. 


B. Survival probability distributions and average transition 
times in finite size conditions 

As described in section II, sets of 50 to 100 simula¬ 
tions were carried out at five different supersaturation 
levels S 1 -S 5 (see I). From each set of simulations an em¬ 
pirical survival probability distribution (ESP) has been 
constructed. The average transition time in finite-size 
conditions tn has been computed for each supersatura¬ 
tion level by a non-linear least square fitting of the ESP 
with the expression Pq = exp (—t/rjv), hereafter referred 
to as the theoretical survival probability (TSP). Evalu¬ 
ating the statistical compatibility between the ESP and 
the TSP with the protocol described in Ref. 25, allowed 
ensuring that the crucial hypothesis of negligible bias de¬ 
position at the transition state has been satisfactorily 
fulfilled for all the simulation sets S 1 -S 5 . In Fig. 6 a both 
the ESP constructed from WTmetaD simulations and the 
fitted TSP are reported; in Fig. 6 b the supersaturation- 
dependent values of the finite-size transition time tn is 
reported as a function of supersaturation, together with 
the average acceleration factor a (in the inset). Fig. 6 b 
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FIG. 6 . a) Survival probability distributions obtained for different supersaturations (*Si=11.4, £ 2 =8.7, 83 = 6 . 8 , 84 = 6 . 0 , 55 = 5 . 4 ). 
b) Average nucleation times in finite size systems tn■ The average acceleration factor a as a function of supersaturation 
is displayed as an inset, c) Nucleation rates calculated from finite-sized WTmetaD simulations (Jn) and rescaled to the 
macroscopic limit J. The blue dashed line represents the result of the fitting of the Jn data with Eq. 14. Computed values of 
J, Jn, tn, ol, convergence of the tn values and their associated errorbars are reported in the Supplementary Information (SI). 


highlights how the application of WTmetaD allows to di¬ 
rectly simulate nucleation events characterised by transi¬ 
tion times of the order of 10 4 seconds, thus significantly 
expanding the range of transition times that could be 
reached with standard MD in a similar simulation setup. 

C. Nucleation rates 

As discussed in section IID nucleation rates in the 
confined regime Jn can be directly computed as Jn = 
(t n V)-\ In Fig. 6 c values of Jn are reported as a 
function of the supersaturation 5, clearly showing a sig¬ 
nificant extension of the accessible nucleation rates do¬ 
main, which is increased up to ten orders of magnitude 
compared to the domain typically accessible by standard 
molecular dynamics 12,23,24,33,34 . 

The surface tension 7 obtained from the fitting of 
the Jn values computed from WTmetaD corresponds to 
7 = 16.9 mN/m, a value that nicely extrapolates the 
data of Goujon et a /. 35 for the same system in the tem¬ 
perature range between 85 K and 135 K as shown in the 
Supplementary Information (SI). The surface tension has 


been considered independent from 5, as typically done in 
CNT. We have found that this choice allows to well de¬ 
scribe the Jn data directly computed from simulations 
while keeping at a minimum the number of fitting param¬ 
eters. We have also verified that considering 7 a linear 
function of supersaturation does not noticeably improve 
the description of the WTmetaD data. 

In Fig. 7, it can be seen that the finite size correction 
is negligible at high supersaturation (S > 11.4), where 
the nucleation barriers AF* and Aare almost in¬ 
distinguishable. However, for lower supersaturation lev¬ 
els reaches values accounting for up to two orders of 
magnitude of difference between J and Jn • This finding 
highlights the importance of properly handling finite size 
effects when investigating transitions at low supersatu¬ 
rations. A breakdown of the contributions of the terms 
appearing in (j) is also reported in Fig. 7. It can be seen 
that, as expected, the contribution of the term /*//^ 
negligible over the entire supersaturation domain. De¬ 
spite the term Z/Zn having a slightly heavier impact on 
0 , it can be seen that the finite-size correction is substan¬ 
tially captured by considering only the exponential term 
in Eq. 15. 
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FIG. 7. Breakdown of the contributions to the finite size 
correction 0 of the factors exp(/3(AF^ — AF*)), ///w, and 
Z/Z*. 

Nucleation rates in macroscopic conditions J are thus 
computed as: 

J = (t>J N ~ J N exp (0 (AF£ - AF*)) (20) 

and reported in Fig. 6c. Nucleation rates rescaled ex¬ 
plicitly accounting also for the term Z/Zn are reported 
in the SI. 


VI. CONCLUSION 

To conclude, in this work we have shown that WT- 
metaD can be applied to the direct calculation of nucle¬ 
ation rates, proving to be particularly useful to tackle 
the timescale limitations that plague small-scale nucle¬ 
ation simulations. In the case of argon condensation this 
implies being able to simulate nucleation in fairly small 
systems, while efficiently reaching timescales of the or¬ 
der of 1 x 10 4 s with ordinary computational resources. 
Moreover, we have highlighted that rate calculations from 
small scale simulations require a systematic assessment 
of finite-size effects. Being able to simultaneously ad¬ 
dress both timescale and finite-size limitations allows to 
significantly extend the range of nucleation conditions 
that can be directly investigated with computationally 
efficient, small-scale molecular simulations. 
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I. TRANSITION TIMES: TABLES AND CONVERGENCE 

Average transition times tat, rates in the confined system Tat, and corrected, macroscopic rates rates J are reported 
in Tab. SI. In order to measure the compatibility between the empirical survival probability distribution computed 
from WTmetaD simulations and the exponential distribution theoretically predicted by the law of rare events a 
Kolmogorov-Smirnov (KS) test has been carried out 1 . The p-value associated to the KS statistic is reported in 
Tab. SI. The p-value, in this context, represents the probability that the survival probability distribution obtained 
from WTmetaD obeys the theoretical exponential distribution. In all cases the p-value is well above the standard 
significance level of 0.05. In Ref. 1 this analysis approach is validated and discussed at length. 


Label 

Tn 

[s] 

Jn 

[cm- 3 s” 1 ] 

J 

[cm” 3 s” 1 ] 

p-value 

a 

Si 

5.75 ±0.65xl0“ 8 

1.5±0.34xl0 25 

3.04±0.70xl0 25 

0.84 

2.8 

s 2 

2.33 ±0.33x10“® 

2.8±0.82xl0 22 

8.64±2.53xl0 22 

0.67 

1.8xl0 2 

s 3 

8.02 ±1.96xl0“ 2 

6.4±3.3xl0 18 

5.09±2.65xl0 19 

0.38 

2.4x10 s 

s 4 

3.61 ±0.76xl0 1 

1.26±0.56xl0 16 

2.57±1.14xl0 17 

0.56 

6.3xl0 7 

s 5 

3.13 ±0.84xl0 4 

1.30±0.75xl0 13 

1.35±0.78xl0 15 

0.37 

1.7xlO n 


TABLE SI. Nucleation rates extracted from WTmetaD simulation, p-value associated to the Kolmogorov-Smirnov statistic 1 , 
and average acceleration factor. 
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The characteristic time convergence with the number of independent samples been reported in Fig. SI. 


10 

r _, 8 

CQ 


b- 


4 

2 

0.2 
0.15 
" 0.1 
0.05 


Sj 

4 

' CQ ' 3 



55 

h 2 

■i 





40 60 80 

samples 


20 30 40 

samples 

_ x 10 4 


30 40 

samples 


40 60 80 

samples 


s 3 

60 

s 4 


^40 


' ~~ -- 

1 — 1 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

20 



0 




°10 20 30 40 50 

samples 


FIG. SI. Convergence of the characteristic time with the number of independent samples. 


Errorbars associated to the characteristic time are computed as the standard deviation of r obtained from a 
bootstrap-like analysis. The convergence of the errorbars as a function of the number of the bootstrap samples is 
displayed in Fig. S2. 



samples 


FIG. S2. Convergence of the error estimation as a function of the number of bootstrap samples. 
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II. IMPORTANCE OF THE ZELDOVICH FACTORS RATIO 
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FIG. S3. Nucleation rates obtained from WTmetaD simulation ( Jn ), corrected for finite-size effects using Eq. 20 (J), and 
corrected for finite size explicitly accounting for the ratio of Zeldovich factors (Jz). It can be seen that Jz values are well 
within the errorbar of the J values. Literature data 2 obtained with different post-processing techniques at high supersaturation 
has been reported (SP survival probability, MFPT mean first passage time, YM Yasuoka-Matsumoto, DO direct observation) 
together with an independent calculation of the nucleation rate at S=11.4 performed applying the SP method to standard MD 
simulations. 
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III. SURFACE TENSION 


Fitting the WTmetaD characteristic times allows to compute surface tension 7 . In our case the surface tension 
is 7 = 16.9 mN/m. Our estimate agrees well with the trend recently reported in the literature 1 2 3 for a truncated 
Lennard-Jones potential. In order to appreciate such agreement, in Fig S4 the a value computed at 80.7 K is reported 
together with surface tension values of Ref. 3. 
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FIG. S4. Surface tension as a function of temperature. The value calculated from the fitting of the nucleation rates J is 
reported together with the data from Goujon et al. 3 . 
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